Pulsar timing arrays as imaging gravitational wave telescopes: 
angular resolution and source (de) confusion 
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Pulsar timing arrays (PTAs) will be sensitive to a finite number of gravitational wave (GW) 
"point" sources (e.g. supermassive black hole binaries). N quiet pulsars with accurately known 
distances d pu i sar can characterize up to 2N/7 distant chirping sources per frequency bin Af gm = 1/T, 
and localize them with "diffraction limited" precision 88 > (l/SNR)(A Btu /d pu i Bar ). Even if the pulsar 
distances are poorly known, a PTA with F frequency bins can still characterize up to (2JV/7)(1— Jp) 
sources per bin, and the quasi-singular pattern of timing residuals in the vicinity of a GW source 
still allows the source to be localized quasi-topologically within roughly the smallest quadrilateral of 
quiet pulsars that encircles it on the sky, down to a limiting resolution SO > (1/SNR) \J X gw / 'd pu i sar . 
PTAs may be unconfused, even at the lowest frequencies, with matched filtering always appropriate. 



Our Local Group of galaxies is sprinkled with milli- 
second pulsars - natural clocks of extraordinary stability. 
Gravitational waves (GWs) passing through the Milky 
Way, after being generated e.g. by the inspiral of two 
supermassive black holes in a distant galaxy, generate 
fluctuations in the time of arrival (TOA) of the pulses 
at the Earth [HQ. In the future, we are likely to de- 
tect such GWs via their coherent imprint on the TOA 
fluctuations from a collection of pulsars distributed on 
the sky: a "pulsar timing array" (PTA). Much research 
has focused on using PTAs to study stochastic GW back- 
grounds (|3j-|5| and references therein). Recently, various 
authors have begun to study the ability of PTAs to detect 
and characterize individual GW point sources 

Continuing in this direction, this paper is concerned 
with conceptually clarifying the theoretical behavior and 
capabilities of PTAs as GW point source telescopes. We 
address two related issues, (i) A PTA may be sensitive to 
so many GW sources that it becomes "confused" - i.e. 
unable to disentangle and individually characterize the 
sources. When does a PTA become "confusion limited" 
rather than sensitivity limited? How many GW sources is 
it capable of individually characterizing? (ii) When a set 
of GW point sources can be individually characterized, 
how well can their angular positions be determined? 

Regarding issue (i) we will see that PTAs with many 
pulsars can characterize many GW sources per frequency 
bin; the traditional rule of thumb that a gravitational 
wave detector becomes confused when there is more 
than about one source per bin is too pessimistic for 
PTAs. Regarding issue (ii) we must distinguish pulsars 
whose distances are known accurately or poorly relative to 
A gjil /(l+cos#), where X gw is the gravitational wavelength 
and 9 is the angle between pulsar and source. Pulsars 
with accurately known distances can angularly localize 
a GW source very precisely; each such pulsar acts like 
a single baseline of a diffraction-limited radio interfer- 
ometer array - with the radio wavelength replaced by 
the gravitational wavelength, and the length of the radio 



baseline replaced by the distance from the pulsar to the 
Earth! The contribution from pulsars with poorly known 
distances is more interesting: due to a quasi-singularity 
in the pattern of timing residuals near the location of the 
GW source, the source can still be localized surprisingly 
well, for reasons that have less to do with diffraction, and 
more to do with topology! 

Basic Formalism. We will label the 3 spatial direc- 
tions with the latin indices fc,i,m = 1,2,3}, raised 
and lowered with S lJ and Sij. The N pulsars in the net- 
work are labelled by the greek indices {a, j3 = 1, . . . , N}, 
raised and lowered with S a " and 8 a p. We follow the 
Einstein summation convention: repeated indices (one 
upper, one lower) are summed. 

A gravitational wave on Minkowski space is described 
in trans verse-traceless (TT) gauge [3] by the line ele- 
ment ds 2 = —dt 2 + [5^ + 2hij\dx % dx 3 '. In this gauge, 
the x = constant worldlines are timelike geodesies; along 
such worldlines, the proper time r is the coordinate time 
t. To avoid notational clutter, let us start with just a sin- 
gle gravitational plane wave travelling in the h direction: 



hi 3 (t, x) 



rf/^(/)e 2 " /(jV£ - t) ; 



(1) 



it is straightforward to extend the following analysis to 
a sum of to = 1, . . . , M plane waves, each travelling in a 
different direction n m ; this extension is discussed below. 
Throughout this paper, we use "dot product" notation 
to mean contraction with the unperturbed 3-metric 8^: 
a ■ b = 8ija l b 3 ; and hats denote unit 3-vectors: a ■ a = 1. 

If an electromagnetic flash is emitted from position Xi 
at time tj, what is its arrival time t at position Xfl If 
we define Xfi = Xf — Xi = XfiXfi then, at zeroth order 
(i.e. in the absence of gravitational waves) the answer is 
to = U + Xfi. Solving the geodesic equation to first order 
in hij yields the perturbed result t = to + St where: 

r ih ii (f)xlxi\e 2 ™K A - s f- to *>-e 2 ™ f ( il - s '-- t ^} 

St = / df lAJ> fl fl \ i. (2) 

J 2wf{l-n ■ Xfi) 
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Now consider an observer at fixed spatial position x = 
receiving signals from a = 1, ...,JV pulsars at spatial 
positions r a = r a r a . For pulsar a, the TOA fluctuation 
3ta(to), as a function of the unperturbed TOA to, is 



where 



St a (t ) 



st a (f) 



df8t a (f)e 



-2mft 



ihMfirUl-VM)) 



(3) 



(4) 



27r/(l+n-f a ) 

and, for later convenience, we have defined the phase 

V a {f)=e 27rifra{1+fl ' fa) . (5) 

The measured TOA fluctuations s a (to) from pulsar a 
are gravitational wave signal St a (to) plus noise n a (to): 



s a (t ) = 5t a (t ) + n a (t ). 



(0) 



We take the noise to be stationary and gaussian, so it 
is characterized by its correlation function C a p{T) or, 
equivalently, its spectral density S a p(f) — C a p{f): 



C a0 {T) = n a {to + T)n p (t ) (7a) 
S(f - f')S aP (f) = K(f)Mf')- (7b) 

We also take the noise to be uncorrelated between dif- 
ferent pulsars: S a /3(f) = S a (f)5 a fs. Let us define the 
natural noise- weighted inner product: 

/>oo 

(3 (1) l5 (2) H / df~g^{fy[S-\f)]^~gf\f). (8) 



Then matched filtering will detect a given gravitational 
wave signal with expected signal-to-noise ratio squared 
(SNR 2 ) given by 



N 



SNR 2 = (St\6t) = SNR a 



SNR,, 



df 



[SUM 
s a (f) 



(9a) 
(9b) 



where St a (f) is given by (j4|). When a gravitational wave 
signal (which depends on various parameters £ fc ) is de- 
tected with sufficient SNR, the likelihood function may 
be approximated as a gaussian oc cxp[— (l/2)£*Tfcz£ ! ] 
near its peak, and the expected inverse covariance matrix 
is the Fisher information matrix, given by 



r -f dt 



dt_ 



(io) 



We are interested, in particular, in the angular resolution 
of a PTA. Define an orthonormal triad from n and two 



other unit vectors rh^ (fi = 1, 2); let 7 M be the rotation 
angle around . The 2x2 angular part of Tki is 



fd[5t] d[6t] 



-pa 



\d~ff 1 
df 



d[St a (f)} 



N 



1 d[8t a (J)] 



(Ha) 
(Hb) 



To evaluate these angular derivatives, we act with the 



infinitessimal rotation R; 



hj ~ "ij - tijkm'^j' 1 on the grav- 
itational wave field, but not on the pulsar positions: e.g. 
d(l+h-r a )/d~f p ' = tijktif^mh and d{h ij (f)f i a fi]/d-fi 1 = 



^h^^e^r^f^m^. In this way, we find 



d[St a (f)} _ i[A(f)+B(f)) 



2tt/ 



(12) 



where A comes from differentiating the phase V a in 
and B comes from differentiating everything else: 



A = 2irifr a 



B = 



l+n ■ f a 



1 + h ■ r, 
2c' 



f J a n m e k i m 



l+n ■ f 



[l-V a ]. 



(13a) 
(13b) 



In understanding the meaning of these equations, we 
should distinguish two cases: (i) pulsars whose distances 
r a are known accurately relative to X gw / (l + ^'?a)> so V a 
is known; and (ii) pulsars whose distances r a are known 
poorly relative to X gw /(1 + n ■ f a ), so V a is essentially a 
random phase. We consider these two cases in turn. 

Pulsars with accurately known distances. First 
consider a monochromatic gravitational plane wave of 
frequency f gw — c/\ gwi and a pulsar whose distance r a 
is known accurately relative to X gw /(1 + n ■ f a ). Then, 
if 2-irr a /X gw 3> 1, the A term dominates the B term in 
Eq. (fT2")l and we have 



SNR 2 



(14) 



Since the second fraction on the right-hand side of this 
equation is typically 0(1), this says that when a pulsar 
at a well known distance r a registers a gravitational wave 
with signal-to- noise level SNR a , its contribution to T~ 
is typically ~ (27rr Q /A gw ) 2 SNR 2 . In other words, 
each such pulsar acts is just like one of the baselines of 
a radio interferometer array; but, in this analogy, the 
radio waves are replaced by gravitational waves, and the 
baselines are of galactic length scales and extend in all 
three spatial dimensions - a remarkable instrument! 

Now consider multiple GW sources. At the low 
GW frequencies probed by PTAs (where the expected 
GW point sources are supermassive black hole binaries, 
far from final merger) the frequency of each gravita- 
tional plane wave drifts negligibly over the observation 



3 



timescale T ~ 10 yrs; and over the light travel time from 
the pulsars to the Earth, the frequency drift or "chirp" 
is approximately linear: hij(t,x) = Re{hije~ 2 ' Rlx ( T )}, 



ing only kicks in for very small angular separations 
9 a < yj\ gw / (2nr a ): for such small separations, V a (fgw) 
ceases to be a random phase, and [1 — V a (f gw )} vanishes 



where x(r) = / r + ±/r 2 , r = t - h ■ x , and {Ay, / 0j /} « C- In other words, when X gw /{2irr a ) < 1, Eq. g]) 



are constants. The induced timing residuals for pul- 
sar a are a sum of two peaks in frequency space: an 
"Earth term" at frequency /o, and a "pulsar term" at 
frequency /o — fr a {l + n ■ r a ); if / is large enough (i.e. 
for supermassive black hole binaries of sufficiently high 
mass, sufficiently close to merger) these two peaks may 
lie in separate frequency bins [6(. The number of such 
GW sources that may be individually characterized by 
a PTA may be determined via the following counting 
argument. To fully specify the pattern of timing residu- 
als, we must provide the following information: in every 
GW frequency bin, and for each "Earth term" in that 
bin, we give the associated propagation direction fi, fre- 
quency derivative /, and two complex amplitudes (i.e. an 
the amplitude and phase for both polarization modes), 
for a total of 7 real numbers. On the other hand, since 
the angular dependence of Eq. (Q} contains spherical har- 
monics of arbitrarily high angular momentum order, the 
number of independent measurements collected by the 
PTA is simply 2N per GW frequency bin - namely, the 
measured amplitude and phase of the timing residuals, 
for each pulsar, in each frequency bin To com- 

pletely characterize the individual sources, the indepen- 
dent measurements must outnumber the parameters to 
be determined; that is, the PTA can characterize up to 
an average of 2N/7 chirping GW point sources per GW 
frequency bin. For simplicity, this argument neglects 
"boundary effects" coming from GW sources for which 
the "Earth term" lies within the detectable frequency 
range, while the "pulsar term" does not, or vice versa. If 
we assume that all of the GW sources are monochromatic 
(/ = 0), the maximum number that can be characterized 
improves only slightly to 2N/6 per frequency bin, but the 
fitting procedure becomes much easier since we can treat 
each GW frequency bin independently. If the PTA can 
disentangle and characterize the individual sources, one 
expects the angular resolution to be diffraction limited 
59 > {1 / SNR) (\ gw /r p U i sar ) [the more precise expecta- 
tion is given by Eq. (|14|)] . 

Pulsars with poorly known distances. If the pul- 
sar distance r a is poorly known relative to \ gw / (l+n-r a ), 
then V a becomes a random phase containing essentially 
no information; the A term is washed out, and only the 
B term remains in Eq. (JTSJ) . Such pulsars no longer con- 
tribute diffraction-limited information, but all is not lost! 

Let us start by giving the key idea, roughly. In spher- 
ical coordinates, with the GW source at the north pole 
and pulsar a at (0 a ,fa)) the factor hij(f)f a f 3 a / (l+h-f a ) 
in Eq. (jH) is the familiar [HI pattern cx cos2ip a (\ + 
cos# a ), which is singular (since the 9 a — > limit de- 
pends on (fa). Ultimately this singularity is smoothed 
out by the [1 — V a ] factor in Eq. but this smooth- 



is guasi-singular at 8 a =0; it becomes genuinely singu- 
lar in the limit A gw /(2irr a ) — > 0. The fact that 5t a (f) 
varies rapidly (with cos 2ip a dependence) around a tiny 
circle of radius y/ X gw / (2irr a ) < 9 a <C 1 surrounding the 
quasi-singularity is the key to localizing the GW source. 

To see this in more detail, split the pulsar directions f a 
into components parallel and perpendicular to the GW 
source direction z: f a — p a sm 9 a + z cos 9 a . Now ap- 
proach the quasi-singularity in two steps: first consider 
the "weaker" limit y / X gw /(2Trr a ) < 9 a <C 1 in which 9 a 
is small but [1— V a ] is not; then proceed to the "stronger" 
limit 9 a -C \J X gw / (2nr a ) <C 1 in which 9 a and [1 — V a \ 
are both small. In the weaker limit, Eq. (|4]) becomes 



st a (f) « -7M/)#ft[i - *>«(/)] (is) 

"7 



while Eq. (fT^j) becomes 



d[St a (f)} _ 2» C ap (/) 
Sj» nf 9 a 



[1-VM)] (16) 



where 



C an {f) = h l3 (f)p a [eVH + P^n<m™e fe;m ] . (17) 



Thus (still in the weaker limit) we have: 



4 C 'afj,(J 'gw)C 'ais(f ' gw) 



SNR 2 91 IhM^pipLl 2 



(18) 



Since the second fraction on the right-hand side of this 
expression is generically 0(1), this says that when a pul- 
sar is near (but not too near) a GW source on the sky, 
its contribution to is typically - (4/6' 2 l )SNR 2 l . 

In the stronger limit 9 a <C y // X gw /(2nr a ) <C 1, Eqs. 
O and CEU) imply that St a (f) and d[5t a {f)]/ are 
smooth and vanishing at 9 a = 0. So as 9 a decreases, 
initially increases as l/9 a , and then drops to zero; 
in between it attains a maximum value: 



X, 



(19) 



gw 



at a separation angle 9 a m y/X gw / {2-nr a ) 1. 

Now consider multiple GW sources. We can repeat 
the previous section's counting argument, except that 
we must now include the N unknown pulsar distances 
when we are counting the parameters needed to spec- 
ify the pattern of timing residuals; with this modifica- 
tion, we find that a PTA which monitors F different 
GW frequency bins can completely characterize up to 
an average of (2A r /7)(l — 1/2-F) "chirping" sources [or 
(2A/6)(1 — 1/2F) monochromatic sources] per bin. Note 
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FIG. 1: The 4 circles represent 4 pulsars that form a small 
square on the sky. The timing residuals of all 4 pulsars are os- 
cillating with the same amplitude and period T gw , but differ- 
ent phases; each pulsar's oscillation is 180° out of phase with 
its 2 nearest neighbors. (This figure depicts this by showing 
4 different moments in the oscillation cycle: black, white, or 
grey circles indicate that, at that moment, the pulses are ar- 
riving early, late, or "on time," respectively.) This signature 
indicates that the square contains a GW point source. 



that, although we didn't know the pulsar distances a pri- 
ori, they are determined, in principle, by the fit fl7| jljj . 

If the PTA can disentangle and characterize the in- 
dividual sources, how well can they be angularly local- 
ized? To answer this question, one should ask, for each 
combination of pulsar and GW source, whether the fit 
to the timing residuals has determined r a accurately or 
poorly relative to X gw /(1 + h ■ f a ); roughly speaking, if 
r a has been determined accurately then we expect the 
pulsar will contribute diffraction limited angular infor- 
mation as described by Eq. (| 14[) : and if r a has been 
determined poorly then we expect the pulsar will con- 
tribute "quasi-singularity limited" angular information 
for that source, as described by Eqs. (fT8|) and fp~9|) . 
Consider the localization of a GW source when all of 
the pulsar's have poorly known distances; as explained 
above, the quasi-singular pattern of timing residuals im- 
plies that the angular localization will be dominated by 
the pulsars that are close to that source on the celestial 
sphere; in particular, it is roughly set by the smallest 
quadrilateral of pulsars that encircles the source on the 
celestial sphere, down to a limiting angular resolution of 
roughly 56 ~ (1/SNR Q ) X gw / d pu i sar [the more precise 
statement is given by Eqs. (fl~8| and (|19))], To understand 



this behavior, consider the example in Fig. [TJ 

Discussion. In the previous sections, we have at- 
tempted to clarify the limits on the capabilities of PTAs 
and, in particular, how these limits depend on factors 
such as the SNR distribution of the GW sources, the 
number and angular distribution of the pulsars relative 
to the GW sources, the distances to the pulsars and the 
precisions of those distances. Our bounds on the angular 
resolution were obtained by Fisher matrix methods; as 
such, these bounds are saturated for GW sources with 
high SNR, and still provide useful guidelines for mod- 
est SNR sources, but will become "loose" for low SNR 
sources. (Note that the relevant SNR here is the total 
SNR of the source in the PTA, which can be high even 
if the SNR per pulsar is not.) Our limits on the number 



of sources that can be individually characterized relied 
on deterministic counting arguments; again, these limits 
will be saturated for high SNR sources, and loose for low 
SNR sources; but even at low SNR, the essential point 
remains that a PTA with many pulsars can distinguish 
many sources per frequency bin, so the traditional rule of 
thumb for confusion (that a gravitational wave detector 
becomes confused when there is > one source per bin) is 
too pessimistic for PTAs. Interesting problems for future 
work include: (i) "tightening" the angular resolution and 
confusion limits at low SNR; (ii) extending this work to 
GW point sources that are near enough that their wave- 
front curvature is significant (llj : (iii) determining the 
circumstances in which pulsar distance determination by 
GW fitting can compete with more traditional methods 
(i.e. VLBI or timing parallax); (iv) clarifying the statis- 
tics of GW sources which are anomalously well charac- 
terized because they are fortuitously located relative to 
one or several pulsars on the sky; (v) quantifying the 
gain from matched filtering (with quasi-singular filters in 
particular) compared to traditional stochastic correlation 
analysis, even when the PTA appears source confused. 
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If Eq. Q only contained spherical harmonics up to an- 
gular momentum order I < l m ax, the number of inde- 
pendent measurements collected by the PTA would be 
limited by (Imax + l) 2 , the total number of spherical har- 
monics up to this order. 

In the special case that the gravitational wave signal is 
due to a single monochromatic gravitational plane wave, 
this fit is highly degenerate, since each r a is only de- 
termined modulo \ gw /(l + n ■ f a ); but, if the GWs are 
significantly non-monochromatic {e.g. "chirping"), signif- 
icantly non-planar [TlT ] or, most importantly, if there are 
two or more GW signals propagating in different direc- 
tions, this degeneracy is broken. 

After completion of this work, we learned that Ref. [l2l ] 
recently emphasized a similar point. 



